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Transition-metal centers are the active sites for many biological and inorganic chemical reac- 
tions. Notwithstanding this central importance, density-functional theory calculations based on 
generalized-gradient approximations often fail to describe energetics, multiplet structures, reaction 
barriers, and geometries around the active sites. We suggest here an alternative approach, derived 
from the Hubbard U correction to solid-state problems, that provides an excellent agreement with 
correlated-electron quantum chemistry calculations in test cases that range from the ground state 
of Fe2 and Fe^" to the addition-elimination of molecular hydrogen on FeO + . The Hubbard U is 
determined with a novel self-consistent procedure based on a linear-response approach. 



Transition metals are central to our understanding of 
many fundamental reactions, as active sites in naturally- 
existing or synthetic molecules that range from metallo- 
porphyrins and oxidoreductases Q to alkene metathesis 
catalysts to light-harvesting photosynthetic complexes 
[3- Despite this relevance, most electronic-structure 
approaches fail to describe consistently or accurately 
transition-metal centers. Examples include neutral and 
charged iron dimers 0], FeO + ^|, Mn(salen) epoxidation 
catalysts Q, or hemeproteins p|. 

In this Letter, we argue that generalized gradient 
approximations (GGA)|gj augmented by a Hubbard U 
term0, already very successful in the solid state 0, 
Hl| . also greatly improve single-site or few-site energies, 
thanks to a more accurate description of self- and intra- 
atomic interactions. Nevertheless, U is not a fitting pa- 
rameter, but an intrinsic response property: as shown 
by Cococcioni and de Gironcoli [l2(, U measures the 
spurious curvature of the GGA energy functional as a 
function of occupations, and GGA+U largely recovers 
the piecewise-linear behavior of the exact ground-state 
energy. U is determined by the difference between the 
screened and bare second derivative of the energy with 
respect to on-site occupations Aj, = ^ \{ [i is the spin- 
orbital, and / the atomic site)[HJ. While in the origi- 
nal derivation U was calculated from the GGA ground 
state, we argue here that U should be consistently ob- 
tained from the GGA+U ground state itself. This be- 
comes especially relevant when GGA and GGA+U differ 
qualitatively (metal vs. insulator in the solid state, dif- 
ferent symmetry in a molecule) . To clarify our approach, 
we first identify in the GGA+U functional the electronic 
terms that have quadratic dependence on the occupa- 
tions: 
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The first term represents the contribution already con- 
tained in the standard GGA functional, modeled here 
as a double-counting term, while the second term is the 
customary "+U" correction. Therefore, U sc / represents 
the effective on-site electron-electron interaction already 
present in the GGA energy functional for the GGA+U 
ground state when U is chosen to be Uj„. Consistency 
is enforced by choosing \J in to be equal to U sc /. The U 
obtained from linear-response E2 (labeled here U ut) is 
also obtained by differentiating Eq. QJwith respect to A^: 
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where m — 1/ J2i( a l) 2 can be interpreted as an effective 
degeneracy of the orbitals whose population is changing 
during the perturbation (to linear order, SXj — a\ (5A^ 
with < = 1 and = £^ afa'j^). Even if in 

principle U sc / depends on U,- n , we find it to be constant 
over a broad interval, as apparent from Fig. ^ U ou t 
is linear in U,-„ for the relevant range of U„ ~ U sc /. 
Thus, from few linear-response calculations for different 
Uj n ground states we are able to extract the U sc / that 
should be used. 

We employ this formulation in the study of the Fe^~ 
and Fe2 dimers and the addition-elimination reaction 
of molecular hydrogen on FeO + : these are paradig- 
matic cases of the challenges for first principles meth- 
ods to accurately reproduce the many low-lying mul- 
tiplet potential energy surfaces associated with transi- 
tion metals. It has been argued that spin density func- 
tional theory can describe the lowest lying; state of a 
given spatial and spin symmetry 0, but difficul- 
ties remain in obtaining accurate multiplet splittings|l5|. 
Our GGA or GGA+U calculations have been per- 
formed with Quantum-ESPRESSO [Hj; coupled clus- 
ter (CCSD(T)) and B3LYP calculations have been per- 
formed with Gaussian03 [18j . 

The iron dimer has been investigated both theoret- 
ically 0, El E3 and experimentally jH El Hi 
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TABLE I: Multiplet splittings (in eV) for Fe^ and Fe 2 at 
several levels of theory, (a) Ref. |2(i|. 



FIG. 1: Linear-response \J ou t calculated from the GGA+Ui„ 
ground state of 7 A„ Fe2, together with the extrapolated U sc /. 
Uo is Uout calculated for \J in —0. 



The experimental photoelectron spectrum of Fe^ be- 
low 2 eV is remarkably simple - there are only two 
prominent peaks, one at 1.0 eV and a second peak 
0.53 eV above it, corresponding to two allowed transi- 
tions to different neutral Fe2 states |23j. A recent multi- 
reference configuration-interaction (MRCI) study|20] has 
assigned the three experimental electronic states involved 
as 8 E~ for Fe^~ and 9 X~ and 7 £~ for Fe2; more recently, 
CCSD(T) has been shown to be in overall agreement [2lj. 
Importantly, these electronic states are consistent with 
the experimental measurements for the anion (funda- 
mental frequency ujq~ 250±20 cm -1 and bond length 
R e =2.10±0.04 A ), and the two neutral Fe2 states, 
which display similar properties (c<Jo=300±15 cm -1 and 
R e =2.02±0.02 A 

We first apply our approach to Fe2 and Fe^~ . We obtain 
a Uo of 2 eV (i.e. when calculated from the GGA ground 
state) and a U sc / of 3 eV (since energies at different U 
are not directly comparable, we average Uo and U sc / over 
all states considered). GGA+U sc f shows a striking and 
consistent agreement with MRCI 20] and our CCSD(T) 
results, correctly identifying both the lowest anion state 
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i g , u.oo cv auuvc. The lowest, singly ion- 
ized neutral states, which differ from Fe^~ only by the 
loss of the spin down or spin up er*(4s) orbital, are 9 £~ 
and 7 Eg . The 9 £~ -> 7 £~ GGA+U sc/ splitting of 0.6 
eV compares very well with theoretical results (MRCI 
and CCSD(T)) and the experimental splitting (0.53 eV) 
in Table The structure of these two states (see Table 
llTjl is also consistent with experimentally observed close 
similarity of R e and ujq for the two neutral states and the 
modest decrease in R e (0.08 A ) and increase in ujq (~ 
50 cm -1 ) with respect to Fe^~ [32{ |. 

In stark contrast with MRCI, CCSD(T) and 
GGA+U SC/ , GGA favors the 8 A g Fe^ state (3d 14 : 
4s 3 : o-qcr*) by as much as 0.9 eV rela- 
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TABLE II: Bond lengths ( A ) and harmonic frequencies, 
u) e , (cm -1 ) for Fe^~ and Fe2, compared to experiment (here, 
fundamental frequencies, loq)- (a) Ref. |2jJ. (b) Ref. 



(3d 13 4s 3 ) which result from the loss of <7*(4s) and cr g (3d) 
electrons, respectively. In addition, these two states have 
differing bond lengths (R e of 1.99 and 2.26 A ) and fre- 
quencies (oj e of 413 cm -1 and 285 cm -1 ), and thus are 
not compatible with experiment 0, [2^] . 

Our second test case explores the potential energy 
surfaces of the highly exothermic (AH < —1.6 eV ) 
addition-elimination reaction of molecular hydrogen on 
bare FeO + . This spin-allowed reaction occurs with ex- 
ceedingly low efficiency (1 in every 100-1000 gas-phase 
collisions results in products), yet when it does proceed 
it is observed to be barrier less [2^, |2|| . This ap- 
parent contradiction has been explained by a two-state- 
reactivity model 0, Eil E(| , wherein the steep reaction 
barriers along the spin surface of the reactants and prod- 
ucts (sextets in both cases) preclude an efficient, exother- 
mic reaction. Instead, the reaction must occur along 
a shallow but excited spin surface (here, the quartet), 
and the reaction bottleneck is the coupling of the two 
surfaces which permits the necessary spin-inversion at 
the entrance and exit channels. For several exchange- 
correlation functional, (including B3LYP1 [jl Ef^ . the re- 
action coordinates have failed to agree qualitatively with 
experiments E5I E riL Et1|. higher level correlated-electron 
calculations |28l l3(ll . o r with the established paradigm of 
a two-state model|29j. 

For the bare FeO + reactant, GGA predicts a 6 S + 
ground state and two nearly degenerate low-lying quartet 



tive to other methods. Neutral states arising from single states, A and 0.84 eV above. GGA+U SC / (5.5 eV) 
ionization of the 8 A g state are 7 A„ (3d 14 4s 2 ) and 9 A g preferentially stabilizes 4 <I> FeO + and yields a 6 E + — > 4 <I> 
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TABLE III: Equilibrium bond lengths, R e (A), harmonic fre- 
quencies, LUe (cm -1 ), and anharmonicities, uj e x e (cm -1 ) for 
the 6 E+ and 4 $ states of FeO+. 




FeO + /H 2 lnt-1 TS-1 lnt-2 TS-2 lnt-3 F e + /H 2 



FIG. 2: Potential energy surface and geometries for the 
FeO + + H2 reaction using GGA (blue) as compared against 
a CCSD(T) reference (black). Solid indicates sextet while 
dashed indicates quartet. (Color online.) 



splitting of 0.54 eV in quantitative agreement with the 
symmetry and splitting (0.57 eV) predicted by CCSD(T). 
The U correction also reduces the 3d character of minor- 
ity spin 7r molecular orbitals which dramatically improves 
bond lengths, harmonic frequencies, and anharmonici- 
ties, as shown in Table ITTTl 

We thus proceed to study the full sextet and quar- 
tet potential energy surfaces (PES) for this reaction. 
We stress that, as is commonly found for open-shell 
transition-metal molecules, several low-lying PES exist 
for each multiplicity and we present results for the lowest- 
lying symmetry of each multiplicity. The U sc / applied in 
this global PES is 5 eV, very close to the average of the 
U sc / (4.93 eV) calculated for the quartet (5.02 eV) and 
sextet (4.84 eV) at each stationary point; the values of 
Uo are similar (quartet = 4.71 eV; sextet = 4.76 eV). 
Although most states possess a U sc / close to the global 
average, the few deviations will be highlighted later. 

Our GGA results for the intermediates (Int) and tran- 
sition states (TS) along the reaction coordinate confirm 
the previously noted failures. Aside from the overesti- 
mate of FcO + splittings, the most notable deviations are 
unusually steep barriers (0.54 eV) along the quartet sur- 
face, lack of spin-crossing near the products, and a a dra- 
matic underestimate in the exothermicity, as depicted in 

Fig. EjajJ. 




FeO + /H 2 lnt-1 TS-1 lnt-2 TS-2 lnt-3 F e + /H 2 



FIG. 3: Potential energy surface and geometries for the FeO + 
+ Eh reaction using GGA+U (5 eV) (blue) as compared 
against a CCSD(T) reference (black). Solid indicates sextet 
while dashed indicates quartet, as in Fig. (Color online.) 
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TABLE IV: Multiplet splittings (in eV) using GGA, GGA+U 
(U = 5 eV except in parentheses, U i n t-z,av= 3.5 eV) and 
CCSD(T). 



With GGA+U (5 eV), we obtain consistency with 
CCSD(T), as shown in Fig. The reactant FeO+ split- 
ting is reduced, the splitting at Int-1 increases, corre- 
sponding to a shallow quartet reaction coordinate, and 
the exothermicity and spin crossover near the products 
are in good agreement with experiment and theoretical 
paradigm0. The quantitative accuracy of GGA+U be- 
comes fully evident in the intermediate splittings (Table 
IIV|I . forward and back reaction barriers (Table [Vj, and 
overall mean absolute errors (MAE) in multiplet split- 
tings that are reduced (with respect to CCSD(T) ref- 
erence) from 0.20 eV for GGA to 0.04 eV for GGA+U. 
Geometries are also improved: the MAE for bond lengths 
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TABLE V: Comparison of GGA, GGA+U (U = 5 eV except 
in parentheses, U4 S = V311 ~ 4 eV) and CCSD(T) forward 
and back reaction barriers (in eV). 
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are reduced from 4.3 pm (GGA) to 2.2 pm (GGA+U). 
0. The GGA+U and CCSD(T) states also possess con- 
sistent orbital occupations and symmetry. 

The few examples of U sc / deviating from 5 eV are pri- 
marily at the exit channel, where large changes in hy- 
bridization occur. The quartet Int-3 is the only case for 
which we obtain a low V sc f (2 eV) which originates from 
the reduced hybridization of Fe 3d states. We chose to 
recalculate the splitting with a U sc / ja ^ that was a lo- 
cal average on the Int-3 states. With this U of 3.5 eV, 
we obtain a splitting of 0.12 eV, in even closer agree- 
ment with CCSD(T). While this reduced hybridization 
of the 3d states is unusual, we stress that it is consis- 
tently predicted in our linear-response approach. Along 
the sextet surface, the iron valence occupations corre- 
spond to 3(i 6 4s 1 , and we find that the the interplay of 3d 
and 4s states to be critical for describing the second bar- 
rier along the sextet reaction surface. A matrix extension 
of our formalism^J considers also the response of the 4s 
orbitals, and we obtain U4 SiSC / = 4.0 eV and U3d jSC /=4.0 
eV around the barrier (U4 S is instead found to be nearly 
zero elsewhere) . Inclusion of the 4s response for both sex- 
tets Int-2 and Int-3 increases the forward reaction barrier 
to 1.16 cV while the backward barrier remains unchanged 
- in accordance with CCSD(T). 

In conclusion, we have shown how a self-consistent 
GGA+U approach can provide a dramatic improvement 
to the description of multiplet potential energy surfaces 
for transition-metal complexes that are otherwise poorly 
described by common exchange-correlation functionals, 
while preserving the very favorable computational costs 
and scaling of local density-based functionals. These 
improvements include spin energetics, state symmetries, 
and quantitative description of complex reaction coordi- 
nates. U has been treated as an intrinsic, non-empirical 
property of the system considered, and never as a fit- 
ting parameter, and it has been obtained through a self- 
consistent extension to the linear-response formulation 
of Cococcioni and de Gironcoli [lij. Such development 
will allow large-scale and accurate calculations |3lJ on 
transition-metal complexes, with applications in the field 
of catalysis, biochemistry, and environmental science. 
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